Run Latent Class Growth Analyses for RAMP & COPING longitudinal analyses of common mental health symptoms: GAD and PHQ, and MASQ (due to CFA indicating tripartite model not a good fit for data
K L Purves
06 August, 2021
This workbook centralises analyses of the PHQ-9, GAD-7 and MASQ total scores longitudinal analyses in the RAMP and COPING samples - identifying the best fit average group trajectories (Latent class growth analyses; LCGA). This corresponds with analytic step 1 in the pre-registration for project “Moderators of symptom trajectories of depression and anxiety during the COVID-19 pandemic”.
The analytic steps that will be followed for PHQ, GAD and MASQ are:
Step 1. Identifying overall group trajectories of anxiety and depression (Latent Growth Curve Modelling; using total scores across all time points for whole sample) Step 2. Assess fit and stability of models Step 3. Identify and conclude the best fitting model Step 4. visualise all relevant trajectories alongside observed data-based trajectories.
Setup
note: some of the install and update processes require command line verification in the console. If you are running this script for the first time and it does not proceed, open and check the console
colour paletttes and plot label list
RAMP colours (first) from website and second is a generated left tetradic (rectangular) colour complement of the RAMP logo green for contrast and complement.
Step 1. Latent Growth Curves
Identify the best fitting latent growth curve form, i.e. that which fits all of the data best, on average, for each outcome.
Notes on timepoints included
The MPlus input time interval specifications go from increasing by 1 for the first 8 follow-ups (we measured with two week intervals) to increasing by 2 for the remaining timepoints when we moved to monthly. The choice of interval is arbitrary, but the relative difference should be retained.
We dropped timepoints 1 & 2 (see ‘Overview of model alterations and adaptations’ below)
Overview of of models
- Linear
- Quadratic
- Piecewise: four pieces, corresponding with lockdowns and easing of restrictions between April 2020 and April 2021 (see below)
- Each of these models will be repeated, including pairwise correlations between consecutive timepoints
Piecewise trajectory key dates
Four pieces
These dates are based on guidelines in England - TO CHECK no major discrepancies across Scotland; NI; Wales
- Lockdown 1 (23 Mar 2020) until first easing of restrictions (23 Jun 2020)
Slope includes data points: Baseline (7 Apr 2020) - Follow-up 5 (16 Jun 2020)
- Summer easing of restrictions (23 Jun 2020) until beginning of Lockdown 2.0 (31 Oct 2020)
Slope includes data points: Follow-up 6 (30 Jun 2020) - Follow-up 11 (20 Oct 2020)
- Winter lockdowns (31 Oct 2020) until Spring 2021 easing of restrictions (8 Mar 2021)
Slope includes data points: Follow-up 12 (17 Nov 2020) - Follow-up 15 (9 Feb 2021)
- Spring 2021 easing of restrictions (8 Mar 2021) until end of study period (6 Apr 2021)
Slope includes data points: Follow-up 16 (9 Mar 2021) - Follow-up 17 (6 Apr 2021)
Notes on lockdowns
England:
Lockdown 2.0: 31 Oct 2020-1 Dec 2020
Easing of restrictions + tier system: 2 Dec 2020 (Tier 3) - 20 Dec 2020 (Tier 4 introduced)
Lockdown 3: 6 Jan 2021 - 8 Mar 2021
Scotland:
Wales:
Northern Ireland:
Covariance coverage issue
After manual checks of each stage, an issue with covariance coverage was detected. The resolution is detailed here.
The covariance coverage (% of participants common to any two time points) between some time points is lower than MPlus default (0.1; 10%). This prevents algorithm for missing data being applied (MLR; maximum likelihood with robust standard errors) and standard errors cannot be computed, even when setting the coverage expectations to zero.
Looking at our data structure and MPlus output, this is due to missing data for follow-up 2 & 3 (RAMP only, no COPING data) resulting in perfect zero covariance coverage with follow-up 17 (currently COPING only; RAMP data extraction issue).
CHANGE TO MODEL: Now running without timepoints 2 & 3
PHQ-9: Current Depression Symptoms Latent Growth Curve Model
Data summary
phq_check <- readRDS("../../../data_clean/phq/phq.clean_merged_total_scores.rds")
# get a list of the total score columns
total_cols <- grep("total_Wave", names(phq_check))
#count how many columns we have
n_cols <- length(total_cols)
# create a list of how many NAs are missing across these 18 columns per row
list_n_missing <- apply(phq_check[total_cols],1,function(x) sum(is.na(x)))
# count how many people are missing data for all columns (NA = total_cols)
list_all_missing <- sum(list_n_missing ==n_cols)
print(dim(phq_check)[1])[1] 52416
[1] 10410
Important note: all numbers correspond exactly with MPlus output numbers (n per wave and n missing all observations etc).
run all models
Don’t rerun models if there is already output for it (to save computational time), do show the output in the console IN FINAL MAKE SURE SCRIPT STATES eval=FALSE so that MPlus models aren’t rerun
Examine model output
This reads in all output, including the gh5 plot files
Check warnings and error messages
- Warnings: check through warnings from console
- Errors: error messages when running models and how each was resolved
Warnings
at the moment this doesnt say which model the warning is from, so any warning requires checking to identify its source. Though order in list will give some indication
[1] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN4 have been fixed to 0", "because the variance of LIN4 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[2] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "2 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10410", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
Warnings are all about missing data (missing on all obs, this number lines up with our expectations from the dataset in the section above)
Otherwise, no warnings. Check one complete.
Note that there are more people missing on all observations having dropped time points 2,3, and 17. Was 7866, but this version is 8217 .
errors
All resolved. See below sections for descriptions.
Check any errors with running any of the models. Output here will not say which model it belongs to, so will need to manually examine if errors occur,
character(0)
Issue 1: Models fail to finish running and converge (resolved)
Some of the models have failed to converge:
- Piecewise with pairwise residual correlations
- Quadratic
- Quadratic with pariwise residual correlations
Steps:
* Increased number of iterations to 5000 (Mplus default is 1000).
* This worked somewhat for the piecewise and quadratic (without pairwise correlations) which now have more informative errors
* This did not work at all for quadratic (with pairwise residual correlation). Increased iterations to 10000 for this one.
* This worked.
Issue 2: Models are no specified due to time point 3 parameter (resolved)
Now all three run, but could not compute standard errors as the models were misidentified or did not converge.
Issues are all with specified parameters now:
quadratic: Problem involving parameter 2 (T_3) * examining the data, I have made an error with the dataset. There are 4325 in this time point. This corresponds to the second follow up that is meant to be dropped from this datafile. Interestingly, T4 onwards is correct. T3 should have 2662 observations. Appears I have wrongly dropped this wave instead of follow up 2. Note: this is exactly what I had done. Corrected and rerun all models
Issue 3: Models are no specified due to low coverage covariacne between t 3 and timepoint 6,8,16 (resolved. CHECK WITh KATIE)
Checked this manually, and there are only 4 people who did both timepoint 3 and timepoint 6, 4 who did timepoint 3 and tomepoint 8, and only 3 who did time point 3 and timepoint 16
- issues to consider here include the fact that having low observed data can caus instability and may lead to ore groups than there really are at a later stage
- On the other hand, the fixed dataset now has good coverage for all other pairings ( min 0.007, representing > 100 people in all cases). Thus dropping this time point might lose more information than it gains.
Options are either dropping time point 3, or dropping coverage requirement.
- tried dropping coverage requirement to see if there remained problems with parameter 3 that would indicate instability.
Issue 4: Quadratic models misspecifed (resolved)
The two quadratic models are not converging. Without pairwise correlations the parameter MPlus flags is the intercept. With pairwise correlation it flags the correlations between time points 15 and 16
- Checked the quadratic output for intercept and noticed some very high variances. I thus specified the model allowing the intercept variance to be freed instead of constrained
- Increasing the number of start value attempts from default 10 2 to 50 5 (50 initial stage starts, 5 final stage optimisations), then 500 initial starts did not work
- Fixing the intercept to 0,40 (variance) 7 (mean) did not work
- made time interval references smaller: MPlus course taught us that models can fail to specify if the numbers that identify relative time intervals are too large. I made it so that 0.5 represented an interval of 2 weeks, instead of 1 THIS WORKED
Check indicators of model success and stability
Tabs to check that:
1. Any residual variance is positive. Significant negative residual variance would indicate a model error. 2. The loadings of indicators and intercepts are 1 3. Estimates late correctly to specific time points (i.e. the relative time references for each time point make sense and are accurate). 4. The intercept is zero for all time points 5. Model estimated correlations between time points look sensible
Output from each relevant file will be printed in each tab for review without needing MPlus software. You may need to scroll across to see full output.
Residuals
check if any residual variances are negative and significant (would indicate model error)
Linear
All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_linear.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Linear with pairwise correlations
All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Quadratic
All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_quadratic.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Quadratic with pairwise correlations
All variances are positive
model <- dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Loadings & Estimates
Check loadings of indicators and intercepts are 1, and estimates relate to specified time points (e.g. lin 0, 1, 2; quad 0, 1, 4)
Note that time references are relative. For most models run here, a difference of 1 refers to a gap of 2 weeks. Thus a difference of 2 refers to a gap of 4 weeks. For Quadratic models, a difference of 0.5 indicates a gap of 2 weeks, and a difference of 1 indicates a gap of 4 weeks. 0 is always the start of a trajectory, before we expect any change to occur. E.g. baseline would be 0 in a linear model.
So linear trajectories (given we are not including follow up time points 1 and 2, which occurred during our 2 week assessment phase) should have a time point reference that looks something like: 0, 3, 4, 5, 6, 7, 8, 9 (two weekly assessment time period ends here, four weeekly assessment period begins), 11, 13, 15, 17, 19, 21, 23, 25
Linear
All is fine
dep.traj.all$phq_latentGrowthCurve_linear.out$output[291:326] %>% str_split("#") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations:
All is fine
dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[295:330] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic.out$output[298:353] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[297:352] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise.out$output[341:432] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[345:436] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Intercepts
check intercept 0 for all time points
Linear
All is fine
dep.traj.all$phq_latentGrowthCurve_linear.out$output[427:444] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations:
All is fine
dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[383:400] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic.out$output[366:383] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[407:423] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise.out$output[458:475] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[507:524] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Correlations
Check that the model estimated correlations between all time points looks sensible.
Linear
All is fine
dep.traj.all$phq_latentGrowthCurve_linear.out$output[653:698] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations:
All is fine
dep.traj.all$phq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[746:792] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic.out$output[709:762] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
dep.traj.all$phq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[792:844] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise.out$output[854:906] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
dep.traj.all$phq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[948:1000] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Summary table of all fit indices
Sorted by BIC, as this is the most informative criteria (as per Muthens)
justSummaries <- do.call("rbind.fill",
sapply(dep.traj.all ,"[", "summaries")) %>%
select(c("Title", "Observations", "Parameters", "AIC", "BIC", "CFI", "TLI", "SRMR", "RMSEA_Estimate")) %>%
arrange(BIC)
justSummaries %>%
kbl() %>%
kable_styling()| Title | Observations | Parameters | AIC | BIC | CFI | TLI | SRMR | RMSEA_Estimate |
|---|---|---|---|---|---|---|---|---|
| Depression symptoms in RAMP and COPING piecewise trajectories With pairwise correlations | 42006 | 46 | 1770170 | 1770568 | 0.995 | 0.994 | 0.009 | 0.015 |
| Depression symptoms in RAMP and COPING Quadratic Growth Curve With pairwise correlations | 42006 | 39 | 1772577 | 1772914 | 0.989 | 0.988 | 0.013 | 0.022 |
| Depression symptoms in RAMP and COPING Linear Growth Curve With pairwise correlations | 42006 | 36 | 1774121 | 1774432 | 0.985 | 0.984 | 0.016 | 0.025 |
| Depression symptoms in RAMP and COPING piecewise trajectory No pairwise correlations | 42006 | 36 | 1774531 | 1774842 | 0.985 | 0.984 | 0.010 | 0.025 |
| Depression symptoms in RAMP and COPING Quadratic Growth Curve No pairwise correlations | 42006 | 25 | 1781344 | 1781560 | 0.968 | 0.970 | 0.017 | 0.035 |
| Depression symptoms in RAMP and COPING Linear Growth Curve No pairwise correlations | 42006 | 21 | 1786100 | 1786282 | 0.956 | 0.960 | 0.021 | 0.040 |
Plots
I show the plots for all runs. Output for all runs is available in the GMM folder on GitHub to review input and output data that generated these plots.
- Observed (dashed lines) and estimated (solid lines) means for best fit trajectory only
- Observed and estimated means: Showing model estimated means from linear, quadratic, and piecewise trajectories (solid lines) overlaid on the observed means from our sample (dashed line)
- Plot estimated mean trajectories for every model run (with and without pairwise correlations)
- plot estimated mean trajectories for the best of each pair of model (i.e. only including te version of the model where contiguous time point residuals were allowed to correlate as these were the best fit for the data)
Best model: Observed and estimated means
Showing model estimated means from piecewise trajectories with pairwise correlations between contiguous time points (solid lines) overlaid on the observed means from our sample (dashed line)
growth_curves_best <- ggplot() +
geom_line(data = dat_dep_best,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_dep_best,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_best
Observed and estimated means
Showing model estimated means from linear, quadratic, and piecewise trajectories (solid lines) overlaid on the observed means from our sample (dashed line) Only showing trajectories that included pairwise correlations between residuals for contiguous time points. This was the best fit for the data.
growth_curves_run2c <- ggplot() +
geom_line(data = dat_dep_run_2c,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_dep_run_2c,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2c
All models
Plotting the estimated means for the trajectories with adjacent residuals correlated as these all had better fit than uncorrelated models
growth_curves_run2 <- ggplot() +
geom_line(data = dat_dep_run_2,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_dep_run_2,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2
Best model from each trajectory type
Show model that performs best of each pair of models (with or without residual correlations between contiguous data points)
Linear: with correlations quadratic: with correlations Piecewise: with correlations
Growth Curve Plot
growth_curves_run2b <- ggplot() +
geom_line(data = dat_dep_run_2b,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_dep_run_2b,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette[c(4, 1:3)]) +
labs(x = "Follow up time point", y ="PHQ9 Score") +
scale_y_continuous(expand = c(0,0), limits = c(6,10), breaks=seq(6, 10, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2b
GAD7: Current Anxiety Symptoms Latent Growth Curve Model
Data summary
gad_check <- readRDS("../../../data_clean/gad/gad.clean_merged_total_scores.rds")
# get a list of the total score columns
total_cols <- grep("total_Wave", names(gad_check))
#count how many columns we have
n_cols <- length(total_cols)
# create a list of how many NAs are missing across these 18 columns per row
list_n_missing <- apply(gad_check[total_cols],1,function(x) sum(is.na(x)))
# count how many people are missing data for all columns (NA = total_cols)
list_all_missing <- sum(list_n_missing ==n_cols)
print(dim(gad_check)[1])[1] 52414
[1] 10968
Examine model output
this reads in all output, including the gh5 plot files
Check warnings and error messages
- warnings
- that none of the residual variances are negative and significant
- loadings of indicators and intercept should all be 1, and estimates should relate to the time points you specified (e.g. lin 0, 1, 2; quad 0, 1, 4)
- intercept 0 for all and correlations look sensible
warnings
at the moment this doesnt say which model the warning is from, so any warning requires checking to identify its source. Though order in list will give some indication
[1] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN4 have been fixed to 0", "because the variance of LIN4 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[2] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "2 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 10969", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
Warnings are all about missing data (missing on all obs, this number lines up with our expectations from the dataset in the section above)
Otherwise, no warnings. Check one complete.
errors
No errors occurred. Used the same alterations as for PHQ (see above)
Check any errors with running any of the models. Output here will not say which model it belongs to, so will need to manually examine if errors occur,
character(0)
Check indicators of model success and stability
Tabs to check that:
1. Any residual variance is positive. Significant negative residual variance would indicate a model error. 2. The loadings of indicators and intercepts are 1 3. Estimates late correctly to specific time points (i.e. the relative time references for each time point make sense and are accurate). 4. The intercept is zero for all time points 5. Model estimated correlations between time points look sensible
Output from each relevant file will be printed in each tab for review without needing MPlus software. You may need to scroll across to see full output.
Residuals
check if any residual variances are negative and significant (would indicate model error) ###### Linear All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_linear.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Linear with pairwise correlations
All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Quadratic
All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_quadratic.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Quadratic with pairwise correlations
All variances are positive
model <- anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Loadings & Estimates
Check loadings of indicators and intercepts are 1, and estimates relate to specified time points (e.g. lin 0, 1, 2; quad 0, 1, 4)
Note that time references are relative. For most models run here, a difference of 1 refers to a gap of 2 weeks. Thus a difference of 2 refers to a gap of 4 weeks. For Quadratic models, a difference of 0.5 indicates a gap of 2 weeks, and a difference of 1 indicates a gap of 4 weeks. 0 is always the start of a trajectory, before we expect any change to occur. E.g. baseline would be 0 in a linear model.
So linear trajectories (given we are not including follow up time points 1 and 2, which occurred during our 2 week assessment phase) should have a time point reference that looks something like: 0, 3, 4, 5, 6, 7, 8, 9 (two weekly assessment time period ends here, four weeekly assessment period begins), 11, 13, 15, 17, 19, 21, 23, 25
Linear
All is fine
anx.traj.all$gad_latentGrowthCurve_linear.out$output[290:327] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[293:330] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic.out$output[298:353] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[297:352] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise.out$output[340:432] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[558:648] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Intercepts
check intercept 0 for all time points
Linear
All is fine
anx.traj.all$gad_latentGrowthCurve_linear.out$output[335:352] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[383:400] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic.out$output[366:383] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[410:427] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise.out$output[458:475] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[718:736] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Correlations
Check that the model estimated correlations between all time points looks sensible.
Linear
All is fine
anx.traj.all$gad_latentGrowthCurve_linear.out$output[653:705] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[746:798] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic.out$output[709:760] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[851:798] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise.out$output[854:906] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
anx.traj.all$gad_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[948:1000] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Summary table of all fit indices
Sorted by BIC, as this is the most informative criteria (as per Muthens)
justSummaries <- do.call("rbind.fill",
sapply(anx.traj.all ,"[", "summaries")) %>%
select(c("Title", "Observations", "Parameters", "AIC", "BIC", "CFI", "TLI", "SRMR", "RMSEA_Estimate")) %>%
arrange(BIC)
justSummaries %>%
kbl() %>%
kable_styling()| Title | Observations | Parameters | AIC | BIC | CFI | TLI | SRMR | RMSEA_Estimate |
|---|---|---|---|---|---|---|---|---|
| Anxiety symptoms in RAMP and COPING piecewise trajectory With pairwise correlations | 41445 | 46 | 1658013 | 1658410 | 0.996 | 0.996 | 0.008 | 0.013 |
| Anxiety symptoms in RAMP and COPING Quadratic Growth Curve With pairwise correlations | 41445 | 40 | 1660647 | 1660992 | 0.989 | 0.988 | 0.013 | 0.021 |
| Anxiety symptoms in RAMP and COPING piecewise trajectory No pairwise correlations | 41445 | 36 | 1661334 | 1661645 | 0.988 | 0.987 | 0.009 | 0.022 |
| Anxiety symptoms in RAMP and COPING Linear Growth Curve With pairwise correlations | 41445 | 36 | 1661947 | 1662258 | 0.986 | 0.985 | 0.016 | 0.024 |
| Anxiety symptoms in RAMP and COPING Quadratic Growth Curve No pairwise correlations | 41445 | 25 | 1667841 | 1668056 | 0.971 | 0.972 | 0.018 | 0.032 |
| Anxiety symptoms in RAMP and COPING Linear Growth Curve No pairwise correlations | 41445 | 21 | 1671934 | 1672116 | 0.960 | 0.963 | 0.021 | 0.037 |
Plots
I show the plots for all runs. Output for all runs is available in the GMM folder on GitHub to review input and output data that generated these plots.
- Observed (dashed lines) and estimated (solid lines) means for best fit trajectory only
- Observed and estimated means: Showing model estimated means from linear, quadratic, and piecewise trajectories (solid lines) overlaid on the observed means from our sample (dashed line)
- Plot estimated mean trajectories for every model run (with and without pairwise correlations)
- plot estimated mean trajectories for the best of each pair of model (i.e. only including te version of the model where contiguous time point residuals were allowed to correlate as these were the best fit for the data)
Best model: Observed and estimated means
Showing model estimated means from piecewise trajectories with pairwise correlations between contiguous time points (solid lines) overlaid on the observed means from our sample (dashed line)
growth_curves_best <- ggplot() +
geom_line(data = dat_anx_best,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anx_best,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_best
#### Observed and estimated means
Plot estimated means from model (solid lines) alongside the observed means (dashed lines) from the data.
Using the best models from each model type pairing (i.e. the version with pairwise correlations of residuals between contiguous time points)
growth_curves_run2c <- ggplot() +
geom_line(data = dat_anx_run_2c,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anx_run_2c,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2c
All models
Plotting the estimated means for the trajectories with adjacent residuals correlated as these all had better fit than uncorrelated models
growth_curves_run2 <- ggplot() +
geom_line(data = dat_anx_run_2,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anx_run_2,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2
Best model from each trajectory type
Only show the best of each pair of models (with or without pairwise residual correlations)
Linear: With correlations Quadratic: With correlations Piecewise: With correlations
growth_curves_run2b <- ggplot() +
geom_line(data = dat_anx_run_2b,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anx_run_2b,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette[c(4, 1:3)]) +
labs(x = "Follow up time point", y ="GAD7 Score") +
scale_y_continuous(expand = c(0,0), limits = c(5,7), breaks=seq(5, 7, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,15.5), breaks=seq(0, 15, by = 1),
labels = timepoint_list) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2b
MASQ: Current anhedonia
Measure has been reverse coded, such that higher scores indicate higher levels of negative affect, to correspond with the direction of the GAD-7 and PHQ 9
Will run these trajectories starting from follow up time 3 as this is the first point where we have both CPOING and RAMP data
Data summary
masq_check <- readRDS("../../../data_clean/masq/masq.clean_merged_total_scores.rds")
# get a list of the total score columns
total_cols <- grep("total_Wave", names(masq_check))
#count how many columns we have
n_cols <- length(total_cols)
# create a list of how many NAs are missing across these 18 columns per row
list_n_missing <- apply(masq_check[total_cols],1,function(x) sum(is.na(x)))
# count how many people are missing data for all columns (NA = total_cols)
list_all_missing <- sum(list_n_missing ==n_cols)
print(dim(masq_check)[1])[1] 37151
[1] 1320
Note: this corresponds to the dataset with all time points, thus, different missing on all data pattern from MPlus run with time points 1,2 and baseline (have more missing all when we exclude the RAMP baseline!)
Examine model output
this reads in all output, including the gh5 plot files
Check warnings and error messages
- warnings
- that none of the residual variances are negative and significant
- loadings of indicators and intercept should all be 1, and estimates should relate to the time points you specified (e.g. lin 0, 1, 2; quad 0, 1, 4)
- intercept 0 for all and correlations look sensible
warnings
at the moment this doesnt say which model the warning is from, so any warning requires checking to identify its source. Though order in list will give some indication
[1] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN1 have been fixed to 0", "because the variance of LIN1 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[2] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("All continuous latent variable covariances involving LIN4 have been fixed to 0", "because the variance of LIN4 is fixed at 0.") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
[3] “c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "3 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS") c("Data set contains cases with missing on all variables.", "These cases were not included in the analysis.", "Number of cases with missing on all variables: 4259", "1 WARNING(S) FOUND IN THE INPUT INSTRUCTIONS")”
Warnings are all about missing data (missing on all obs, this number lines up with our expectations from the dataset in the section above)
Otherwise, no warnings. Check one complete.
errors
Check any errors with running any of the models. Output here will not say which model it belongs to, so will need to manually examine if errors occur,
character(0)
Issue 1: piecewise with correlations (resolved)
piecewise with correlations did not run due to issues with the first linear piece. Likely due to few pieces contributing to this section. Restrained variance to 0 for this piece and it ran successfully.
Check indicators of model success and stability
Tabs to check that:
1. Any residual variance is positive. Significant negative residual variance would indicate a model error. 2. The loadings of indicators and intercepts are 1 3. Estimates late correctly to specific time points (i.e. the relative time references for each time point make sense and are accurate). 4. The intercept is zero for all time points 5. Model estimated correlations between time points look sensible
Output from each relevant file will be printed in each tab for review without needing MPlus software. You may need to scroll across to see full output.
Residuals
check if any residual variances are negative and significant (would indicate model error) ###### Linear All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_linear.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Linear with pairwise correlations
All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Quadratic
All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_quadratic.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Quadratic with pairwise correlations
All variances are positive
model <- anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$parameters$unstandardized
model[model$paramHeader == "Residual.Variances",]Loadings & Estimates
Check loadings of indicators and intercepts are 1, and estimates relate to specified time points (e.g. lin 0, 1, 2; quad 0, 1, 4)
Note that time references are relative. For most models run here, a difference of 1 refers to a gap of 2 weeks. Thus a difference of 2 refers to a gap of 4 weeks. For Quadratic models, a difference of 0.5 indicates a gap of 2 weeks, and a difference of 1 indicates a gap of 4 weeks. 0 is always the start of a trajectory, before we expect any change to occur. E.g. baseline would be 0 in a linear model.
So linear trajectories (given we are not including follow up time points 1 and 2, which occurred during our 2 week assessment phase) should have a time point reference that looks something like: 0, 3, 4, 5, 6, 7, 8, 9 (two weekly assessment time period ends here, four weeekly assessment period begins), 11, 13, 15, 17, 19, 21, 23, 25
Linear
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear.out$output[280:314] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[284:317] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic.out$output[400:450] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[287:337] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise.out$output[490:574] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[340:424] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Intercepts
check intercept 0 for all time points
Linear
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear.out$output[322:338] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[367:383] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic.out$output[463:479] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[392:408] %>%
str_split("/t") %>%
kbl() %>%
kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise.out$output[600:616] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[694:710] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Correlations
Check that the model estimated correlations between all time points looks sensible.
Linear
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear.out$output[600:643] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Linear with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_linear_withpairwiseCorrelations.out$output[687:730] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic.out$output[654:697] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Quadratic with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_quadratic_withPairwiseCorrelations.out$output[736:779] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise.out$output[793:836] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Piecewise with pairwise correlations
All is fine
anhed.traj.all$masq_latentGrowthCurve_piecewise_wPairwise_correlations.out$output[887:930] %>% str_split("/t") %>% kbl() %>% kable_styling(bootstrap_options = "striped")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Summary table of all fit indices
Sorted by BIC, as this is the most informative criteria (as per Muthens)
justSummaries <- do.call("rbind.fill",
sapply(anhed.traj.all ,"[", "summaries")) %>%
select(c("Title", "Observations", "Parameters", "AIC", "BIC", "CFI", "TLI", "SRMR", "RMSEA_Estimate")) %>%
arrange(BIC)
justSummaries %>%
kbl() %>%
kable_styling()| Title | Observations | Parameters | AIC | BIC | CFI | TLI | SRMR | RMSEA_Estimate |
|---|---|---|---|---|---|---|---|---|
| Anhedonia symptoms in RAMP and COPING piecewise trajectory With pairwise correlations | 32892 | 40 | 1795036 | 1795372 | 0.984 | 0.982 | 0.032 | 0.032 |
| Anhedonia symptoms in RAMP and COPING piecewise trajectory No pairwise correlations | 32892 | 35 | 1796887 | 1797181 | 0.978 | 0.977 | 0.026 | 0.037 |
| Anhedonia symptoms in RAMP and COPING Quadratic Growth Curve With pairwise correlations | 32892 | 38 | 1797427 | 1797746 | 0.975 | 0.973 | 0.042 | 0.040 |
| Anhedonia symptoms in RAMP and COPING Linear Growth Curve With pairwise correlations | 32892 | 34 | 1799404 | 1799690 | 0.968 | 0.967 | 0.043 | 0.044 |
| Anhedonia symptoms in RAMP and COPING Quadratic Growth Curve No pairwise correlations | 32892 | 24 | 1802466 | 1802668 | 0.958 | 0.961 | 0.044 | 0.048 |
| Anhedonia symptoms in RAMP and COPING Linear Growth Curve No pairwise correlations | 32892 | 20 | 1806998 | 1807166 | 0.943 | 0.948 | 0.046 | 0.056 |
Plots
I show the plots for all runs. Output for all runs is available in the GMM folder on GitHub to review input and output data that generated these plots.
- Observed (dashed lines) and estimated (solid lines) means for best fit trajectory only
- Observed and estimated means: Showing model estimated means from linear, quadratic, and piecewise trajectories (solid lines) overlaid on the observed means from our sample (dashed line)
- Plot estimated mean trajectories for every model run (with and without pairwise correlations)
- plot estimated mean trajectories for the best of each pair of model (i.e. only including te version of the model where contiguous time point residuals were allowed to correlate as these were the best fit for the data)
Best model: Observed and estimated means
Showing model estimated means from piecewise trajectories with pairwise correlations between contiguous time points (solid lines) overlaid on the observed means from our sample (dashed line)
growth_curves_best <- ggplot() +
geom_line(data = dat_anhed_best,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anhed_best,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_best
#### Observed and estimated means
Plot estimated means from model (solid lines) alongside the observed means (dashed lines) from the data.
Using the best models from each model type pairing (i.e. the version with pairwise correlations of residuals between contiguous time points)
growth_curves_run2c <- ggplot() +
geom_line(data = dat_anhed_run_2c,
aes(x = timepoint, y = score, colour = Model,linetype=Type)) +
geom_point(data = dat_anhed_run_2c,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2c
All models
Plotting the estimated means for the trajectories with adjacent residuals correlated as these all had better fit than uncorrelated models
growth_curves_run2 <- ggplot() +
geom_line(data = dat_anhed_run_2,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anhed_run_2,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2
Best model from each trajectory type
Only show the best of each pair of models (with or without pairwise residual correlations)
Linear: With correlations Quadratic: With correlations Piecewise: With correlations
growth_curves_run2b <- ggplot() +
geom_line(data = dat_anhed_run_2b,
aes(x = timepoint, y = score, colour = Model), size = 1) +
geom_point(data = dat_anhed_run_2b,
aes(x = timepoint, y = score, colour = Model, shape=Model), show.legend = F) +
scale_color_manual(values = pres.palette[c(4, 1:3)]) +
labs(x = "Follow up time point", y ="masq (reverse) Score") +
scale_y_continuous(expand = c(0,0), limits = c(25,29), breaks=seq(25, 29, 1)) +
scale_x_continuous(expand = c(0,0), limits = c(-0.5,14.5), breaks=seq(0, 14, by = 1),
labels = timepoint_list_anh) +
theme(panel.grid.major.y = element_line(size = 0.5,
linetype = 'dashed',
colour = "gray"),
axis.text = element_text(colour="black", size = 14),
axis.title = element_text(colour="black", size = 16),
axis.text.x = element_text(angle = 90,
hjust=1),
legend.key = element_blank(),
legend.text = element_text(colour = "black", size = 11),
panel.background = element_blank())
growth_curves_run2b